Code
pacman::p_load(sp, sf, spNetwork, tmap, classInt, viridis, tidyverse, spatstat, sfdep, raster, maptools, arrow, lubridate)Take-home Exercise 01
1.0 Background
Human mobility, the movement of human beings in space and time, reflects the spatial-temporal characteristics of human behavior. With the advancement Information and Communication Technologies (ICT) especially smart phone, a large volume of data related to human mobility have been collected. By using appropriate GIS analysis methods, these data are potentially useful in supporting smart city planning and management.
In Singapore, one of the important source of data related to human mobility is from Land Transport Authority (LTA) DataMall. Two data sets related to human mobility are provided by the portal, they are: Passenger Volume by Origin Destination Train Stations and Passenger Volume by Origin Destination Bus Stops. One of the limitation of these data sets is that their location are biased to either bus stops or MRT/LRT stations. In 2020, another very interesting human mobility data set called Grab Posisi was released by GRAB, one of the largest shared taxi operator in South-east Asia. There are two data sets been released and one of them is for Singapore.
2.0 Objectives
Geospatial analytics hold tremendous potential to address complex problems facing society. In this study, you are tasked to apply appropriate spatial point patterns analysis methods to discover the geographical and spatio-temporal distribution of Grab hailing services locations in Singapore.
3.0 The Task
The specific tasks of this take-home exercise are as follows:
Using appropriate function of sf and tidyverse, preparing the following geospatial data layer in sf tibble data.frames:
Grab taxi location points either by origins or destinations.
Road layer within Singapore excluding outer islands.
Singapore boundary layer excluding outer islands
Using the extracted data, derive traditional Kernel Density Estimation layers.
Using the extracted data, derive either Network Kernel Density Estimation (NKDE) or Temporal Network Kernel Density Estimation (TNKDE)
Using appropriate tmap functions, display the kernel density layers on openstreetmap of Singapore.
Describe the spatial patterns revealed by the kernel density maps.
4.0 The Data
Apstial data
For the purpose of this assignment, Grab-Posisi of Singapore will be used.
Geospatial data
Road data set from OpenStreetMap of Geofabrik download server. The Malaysia, Singapore, and Brunei coverage should be downloaded.
Master Plan 2019 Subzone Boundary (No Sea) from Data.gov.sg.
5.0. Importing Packages
The R packages used in this project are:
sf, a relatively new R package specially designed to import, manage and process vector-based geospatial data in R.
tidyverse for performing data science tasks such as importing, wrangling and visualising data.
spatstat, which has a wide range of useful functions for point pattern analysis. In this hands-on exercise, it will be used to perform 1st- and 2nd-order spatial point patterns analysis and derive kernel density estimation (KDE) layer.
raster which reads, writes, manipulates, analyses and model of gridded spatial data (i.e. raster). In this hands-on exercise, it will be used to convert image output generate by spatstat into raster format.
maptools which provides a set of tools for manipulating geographic data. In this hands-on exercise, we mainly use it to convert Spatial objects into ppp format of spatstat.
tmap which provides functions for plotting cartographic quality static point patterns maps or interactive maps by using leaflet API.
spNetwork, which provides functions to perform Spatial Point Patterns Analysis such as kernel density estimation (KDE) and K-function on network. It also can be used to build spatial matrices (‘listw’ objects like in ‘spdep’ package) to conduct any kind of traditional spatial analysis with spatial weights based on reticular distances.
rgdal, which provides bindings to the ‘Geospatial’ Data Abstraction Library (GDAL) (>= 1.11.4) and access to projection/transformation operations from the PROJ library. In this exercise, rgdal will be used to import geospatial data in R and store as sp objects.
sp, which provides classes and methods for dealing with spatial data in R. In this exercise, it will be used to manage SpatialPointsDataFrame and SpatiaLinesDataFrame, and for performing projection transformation.
6.0 Importing Data in R environment
6.1 Aspatial Data
6.1.1 Import GrabPosisi Data
Code
df <- read_parquet("data/Apstial/GrabPosisi/part-00000-8bbff892-97d2-4011-9961-703e38972569.c000.snappy.parquet")6.2 Geospatial Data
6.2.1 Import Master Plan 2019 subzone Boundary (No Sea) Data
Code
mpsz_sf <- st_read(dsn= "data/Geospatial/MPSZ-2019", layer="MPSZ-2019")Reading layer `MPSZ-2019' from data source
`C:\j00b00\IS415-GAA\Take-home_EX\Take-home_Ex01\data\Geospatial\MPSZ-2019'
using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS: WGS 846.2.2 Import Road data
Code
roads <- st_read(dsn = "data/Geospatial/malaysia-singapore-brunei-latest-free", layer = "gis_osm_roads_free_1")Reading layer `gis_osm_roads_free_1' from data source
`C:\j00b00\IS415-GAA\Take-home_EX\Take-home_Ex01\data\Geospatial\malaysia-singapore-brunei-latest-free'
using driver `ESRI Shapefile'
Simple feature collection with 1759836 features and 10 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 99.66041 ymin: 0.8021131 xmax: 119.2601 ymax: 7.514393
Geodetic CRS: WGS 847.0 Data Wrangling
Before we proceed with our analysis, we need to preprocess the data to the correct format and clean it.
7.1 Assign the correct CRS
7.1.1 mpsz_sf
Code
mpsz3414 <- st_transform(mpsz_sf, crs = 3414)
st_crs(mpsz3414)Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]7.1.2 roads
Code
roads3414 <- st_transform(roads, crs = 3414)
st_crs(roads3414)Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]7.2 Convert data type of pingtimestamp from character to date-time
7.2.1 df
Code
df$pingtimestamp <- as_datetime(df$pingtimestamp)7.2.2 Extract Grab taxi location points by origins
Code
origin_df <- df %>%
group_by(trj_id) %>%
arrange(desc(pingtimestamp)) %>%
filter(row_number()==1) %>%
mutate(weekday = wday(pingtimestamp, label=TRUE, abbr=TRUE),
start_hr = factor(hour(pingtimestamp)),
day = factor(mday(pingtimestamp)))7.2.3 Creating a simple feature data frame from origin_df
Code
origin_sf <- st_as_sf(origin_df,
coords = c("rawlng", "rawlat"),
crs = 4326) %>%
st_transform(crs = 3414)7.3 Extract Singapore boundary layer excluding outer islands
7.3.1 First, filter out the outer islands
Code
outer_island <- mpsz_sf[mpsz_sf$PLN_AREA_N == "SOUTHERN ISLANDS" | mpsz_sf$PLN_AREA_N == "NORTH-EASTERN ISLANDS" | mpsz_sf$PLN_AREA_N == "WESTERN ISLANDS",]Code
plot(st_geometry(outer_island))
Code
outer_island<- st_transform(outer_island, crs = 3414)7.3.3 Filter out the outer island
Code
main_island <- st_difference(st_union(mpsz3414),st_union(outer_island))Code
plot(st_geometry(main_island))
7.4 Filter out roads from roads3414 that have max speed value above 0.
Code
main_road <- roads3414 %>% filter(maxspeed > 0)7.4.1 Find the roads that are only in Singapore excluding outer islands
Code
sg_roads <- st_intersection(main_road, main_island)7.5 Road layer
Code
tm_shape(main_island) + tm_polygons() + tm_shape(sg_roads) + tm_lines(col="red", siz = 0.03)
7.6 Grab origin layer
Code
tm_shape(origin_sf) + tm_dots(siz = 0.03)
7.7 Converting sf data frames to sp’s Spatial* class
spatstat requires the analytical data in ppp object form. There is no direct way to convert a Spatial* classes into ppp object. We need to convert the Spatial classes* into Spatial object first.
Code
origin <- as_Spatial(origin_sf)
mpsz <- as_Spatial(mpsz3414)
main_island <- as_Spatial(main_island)
sg_roads <- as_Spatial(sg_roads)7.8 Converting Spatial* class into generic sp format
Code
origin_sp <- as(origin, "SpatialPoints")
main_island_sp <- as(main_island, "SpatialPolygons")
sg_roads_sp <- as(sg_roads, "SpatialLines")
origin_spclass : SpatialPoints
features : 28000
extent : 3638.685, 50024.92, 25350.05, 49469.41 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 7.9 Converting generic sp format into spatstat’s ppp format
convert the spatial data into spatstat’s ppp object format.
Code
origin_ppp <- as(origin_sp, "ppp")7.9.1 Plot origin_ppp for visualisation
Code
plot(origin_ppp)
Code
summary(origin_ppp)Planar point pattern: 28000 points
Average intensity 2.502667e-05 points per square unit
Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units
Window: rectangle = [3638.69, 50024.92] x [25350.05, 49469.41] units
(46390 x 24120 units)
Window area = 1118810000 square units7.10 Handling Duplicate Point
check the duplication in a ppp object
Code
any(duplicated(origin_ppp))[1] FALSE7.11 Creating owin object
When analysing spatial point patterns, it is a good practice to confine the analysis with a geographical area like Singapore boundary. In spatstat, an object called owin is specially designed to represent this polygonal region.
Code
main_island_owin <- as(main_island, "owin")7.11.1 Display main_island_owin object
Code
plot(main_island_owin)
Code
summary(main_island_owin)Window: polygonal boundary
37 separate polygons (29 holes)
vertices area relative.area
polygon 1 71 5.63061e+03 8.47e-06
polygon 2 10 1.99717e+02 3.01e-07
polygon 3 12667 6.63014e+08 9.98e-01
polygon 4 (hole) 3 -3.41897e-05 -5.14e-14
polygon 5 (hole) 23 -1.99656e+01 -3.00e-08
polygon 6 (hole) 35 -1.38385e+02 -2.08e-07
polygon 7 (hole) 19 -4.39650e+00 -6.62e-09
polygon 8 (hole) 270 -1.21455e+03 -1.83e-06
polygon 9 (hole) 3 -4.95057e-02 -7.45e-11
polygon 10 (hole) 3 -3.65499e-03 -5.50e-12
polygon 11 (hole) 38 -7.79904e+03 -1.17e-05
polygon 12 (hole) 3 -5.99535e-04 -9.02e-13
polygon 13 (hole) 3 -3.04561e-04 -4.58e-13
polygon 14 (hole) 3 -7.43616e-06 -1.12e-14
polygon 15 (hole) 6 -8.37554e-01 -1.26e-09
polygon 16 (hole) 4 -2.86396e-01 -4.31e-10
polygon 17 (hole) 3 -1.81439e-04 -2.73e-13
polygon 18 (hole) 3 -8.68789e-04 -1.31e-12
polygon 19 (hole) 3 -4.46076e-04 -6.71e-13
polygon 20 (hole) 3 -3.39794e-04 -5.11e-13
polygon 21 (hole) 317 -5.11280e+04 -7.69e-05
polygon 22 (hole) 5 -2.92235e-04 -4.40e-13
polygon 23 (hole) 3 -4.52043e-05 -6.80e-14
polygon 24 (hole) 3 -3.90173e-05 -5.87e-14
polygon 25 (hole) 5 -2.44411e-04 -3.68e-13
polygon 26 (hole) 4 -2.18616e-04 -3.29e-13
polygon 27 (hole) 4 -4.28453e-01 -6.45e-10
polygon 28 (hole) 4 -2.54488e-04 -3.83e-13
polygon 29 (hole) 3 -9.59850e-05 -1.44e-13
polygon 30 (hole) 41 -4.01660e+04 -6.04e-05
polygon 31 (hole) 3 -4.14099e-04 -6.23e-13
polygon 32 (hole) 5 -3.64686e-02 -5.49e-11
polygon 33 30 2.80002e+04 4.21e-05
polygon 34 27 1.50315e+04 2.26e-05
polygon 35 285 1.61128e+06 2.42e-03
polygon 36 91 1.49663e+04 2.25e-05
polygon 37 71 8.18750e+03 1.23e-05
enclosing rectangle: [2667.54, 55941.94] x [21448.47, 50256.33] units
(53270 x 28810 units)
Window area = 664597000 square units
Fraction of frame area: 0.4337.12 Combining point events object and owin object
Extract Grab origin point that are in main island in Singapore
Code
originSG_ppp = origin_ppp[main_island_owin]Code
summary(originSG_ppp)Planar point pattern: 27821 points
Average intensity 4.186147e-05 points per square unit
Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units
Window: polygonal boundary
37 separate polygons (29 holes)
vertices area relative.area
polygon 1 71 5.63061e+03 8.47e-06
polygon 2 10 1.99717e+02 3.01e-07
polygon 3 12667 6.63014e+08 9.98e-01
polygon 4 (hole) 3 -3.41897e-05 -5.14e-14
polygon 5 (hole) 23 -1.99656e+01 -3.00e-08
polygon 6 (hole) 35 -1.38385e+02 -2.08e-07
polygon 7 (hole) 19 -4.39650e+00 -6.62e-09
polygon 8 (hole) 270 -1.21455e+03 -1.83e-06
polygon 9 (hole) 3 -4.95057e-02 -7.45e-11
polygon 10 (hole) 3 -3.65499e-03 -5.50e-12
polygon 11 (hole) 38 -7.79904e+03 -1.17e-05
polygon 12 (hole) 3 -5.99535e-04 -9.02e-13
polygon 13 (hole) 3 -3.04561e-04 -4.58e-13
polygon 14 (hole) 3 -7.43616e-06 -1.12e-14
polygon 15 (hole) 6 -8.37554e-01 -1.26e-09
polygon 16 (hole) 4 -2.86396e-01 -4.31e-10
polygon 17 (hole) 3 -1.81439e-04 -2.73e-13
polygon 18 (hole) 3 -8.68789e-04 -1.31e-12
polygon 19 (hole) 3 -4.46076e-04 -6.71e-13
polygon 20 (hole) 3 -3.39794e-04 -5.11e-13
polygon 21 (hole) 317 -5.11280e+04 -7.69e-05
polygon 22 (hole) 5 -2.92235e-04 -4.40e-13
polygon 23 (hole) 3 -4.52043e-05 -6.80e-14
polygon 24 (hole) 3 -3.90173e-05 -5.87e-14
polygon 25 (hole) 5 -2.44411e-04 -3.68e-13
polygon 26 (hole) 4 -2.18616e-04 -3.29e-13
polygon 27 (hole) 4 -4.28453e-01 -6.45e-10
polygon 28 (hole) 4 -2.54488e-04 -3.83e-13
polygon 29 (hole) 3 -9.59850e-05 -1.44e-13
polygon 30 (hole) 41 -4.01660e+04 -6.04e-05
polygon 31 (hole) 3 -4.14099e-04 -6.23e-13
polygon 32 (hole) 5 -3.64686e-02 -5.49e-11
polygon 33 30 2.80002e+04 4.21e-05
polygon 34 27 1.50315e+04 2.26e-05
polygon 35 285 1.61128e+06 2.42e-03
polygon 36 91 1.49663e+04 2.25e-05
polygon 37 71 8.18750e+03 1.23e-05
enclosing rectangle: [2667.54, 55941.94] x [21448.47, 50256.33] units
(53270 x 28810 units)
Window area = 664597000 square units
Fraction of frame area: 0.4337.12.1 Plot originSg_ppp
Code
plot(originSG_ppp)
8.0 First-order Spatial Point Patterns Analysis
Deriving kernel density estimation (KDE) layer for visualising and exploring the intensity of point processes
8.1 Kernel Density Estimation
8.1.1 Compute kernel density estimation using automatic origin points in main land in Singapore
Code
kde_originSG_bw = density(originSG_ppp,
sigma = bw.diggle,
edge = TRUE,
kernel="gaussian")8.1.2 Display the kernel density derived
Code
plot(kde_originSG_bw)
8.1.3 Rescalling KDE Values
Code
originSG_ppp.km <-rescale(originSG_ppp, 1000, "km")Re-run density() using the resale data set and plot the output kde map.
Code
kde_originSG.bw <- density(originSG_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")8.1.4 Display the kernel density derived
Code
plot(kde_originSG.bw)
8.2 Converting KDE output into grid object Convert it so that it is suitable for mapping purposes
Code
gridded_kde_originSG_bw <- as.SpatialGridDataFrame.im(kde_originSG.bw)
spplot(gridded_kde_originSG_bw)
8.2.1 Converting gridded output into raster Convert the gridded kernal density objects into RasterLayer object
Code
kde_originSG_bw_raster <- raster(gridded_kde_originSG_bw)Properties of kde_originSG_bw_raster RasterLayer
Code
kde_originSG_bw_rasterclass : RasterLayer
dimensions : 128, 128, 16384 (nrow, ncol, ncell)
resolution : 0.4162063, 0.2250614 (x, y)
extent : 2.667538, 55.94194, 21.44847, 50.25633 (xmin, xmax, ymin, ymax)
crs : NA
source : memory
names : v
values : -3.709752e-13, 1953.626 (min, max)8.2.2 Assigning projection system in kde_originSG_bw_raster RasterLayer
Code
projection(kde_originSG_bw_raster) <- CRS("+init=EPSG:3414")
kde_originSG_bw_rasterclass : RasterLayer
dimensions : 128, 128, 16384 (nrow, ncol, ncell)
resolution : 0.4162063, 0.2250614 (x, y)
extent : 2.667538, 55.94194, 21.44847, 50.25633 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs
source : memory
names : v
values : -3.709752e-13, 1953.626 (min, max)8.2.3 Visualising the output in tmap
Code
tm_shape(kde_originSG_bw_raster) + tm_raster("v") +
tm_layout(legend.position = c("right", "bottom"),frame = FALSE)
8.3 Comparing Spatial Point Patterns using KDE Compare KDE of childcare at Ponggol, Tampines, Chua Chu Kang and Jurong West planning areas.
8.3.1 Extracting planning area boundary
Code
pg = mpsz[mpsz@data$PLN_AREA_N == "PUNGGOL",]
tm = mpsz[mpsz@data$PLN_AREA_N == "TAMPINES",]
ck = mpsz[mpsz@data$PLN_AREA_N == "CHOA CHU KANG",]
jw = mpsz[mpsz@data$PLN_AREA_N == "JURONG WEST",]8.3.2 Plotting target planning areas
Code
plot(pg, main = "Ponggol")
Code
plot(tm, main = "Tampines")
Code
plot(ck, main = "Choa Chu Kang")
Code
plot(jw, main = "Jurong West")
8.3.3 Converting the spatial point data frame into generic sp format
Code
pg_sp = as(pg, "SpatialPolygons")
tm_sp = as(tm, "SpatialPolygons")
ck_sp = as(ck, "SpatialPolygons")
jw_sp = as(jw, "SpatialPolygons")8.3.4 Creating owin object
Code
pg_owin = as(pg_sp, "owin")
tm_owin = as(tm_sp, "owin")
ck_owin = as(ck_sp, "owin")
jw_owin = as(jw_sp, "owin")8.3.5 Combining origin points and the study area Extract childcare that is within the specific region to do our analysis later on
Code
origin_pg_ppp = originSG_ppp[pg_owin]
origin_tm_ppp = originSG_ppp[tm_owin]
origin_ck_ppp = originSG_ppp[ck_owin]
origin_jw_ppp = originSG_ppp[jw_owin]Next, rescale() function is used to transform the unit of measurement from metre to kilometer.
Code
origin_pg_ppp.km = rescale(origin_pg_ppp, 1000, "km")
origin_tm_ppp.km = rescale(origin_tm_ppp, 1000, "km")
origin_ck_ppp.km = rescale(origin_ck_ppp, 1000, "km")
origin_jw_ppp.km = rescale(origin_jw_ppp, 1000, "km")Plot these four study areas and the locations of the origin points
Code
plot(origin_pg_ppp.km, main="Punggol")
Code
plot(origin_tm_ppp.km, main="Tampines")
Code
plot(origin_ck_ppp.km, main="Choa Chu Kang")
Code
plot(origin_jw_ppp.km, main="Jurong West")
8.3.6 Computing KDE
Code
plot(density(origin_pg_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian"),
main="Punggol")
Code
plot(density(origin_tm_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian"),
main="Tampines")
Code
plot(density(origin_ck_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian"),
main="Choa Chu Kang")
Code
plot(density(origin_jw_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian"),
main="Jurong West")
8.3.7 Storing KDE values
Code
kde_origin_pg.bw = density(origin_pg_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")
kde_origin_tm.bw = density(origin_tm_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")
kde_origin_ck.bw = density(origin_ck_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")
kde_origin_jw.bw = density(origin_jw_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")8.4 Converting KDE output into grid object Convert it so that it is suitable for mapping purposes
Code
gridded_kde_origin_pg_bw <- as.SpatialGridDataFrame.im(kde_origin_pg.bw)
spplot(gridded_kde_origin_pg_bw)
Code
gridded_kde_origin_tm_bw <- as.SpatialGridDataFrame.im(kde_origin_tm.bw)
spplot(gridded_kde_origin_tm_bw)
Code
gridded_kde_origin_ck_bw <- as.SpatialGridDataFrame.im(kde_origin_ck.bw)
spplot(gridded_kde_origin_ck_bw)
Code
gridded_kde_origin_jw_bw <- as.SpatialGridDataFrame.im(kde_origin_jw.bw)
spplot(gridded_kde_origin_jw_bw)
8.4.1 Converting gridded output into raster
Code
kde_origin_pg_bw_raster <- raster(gridded_kde_origin_pg_bw)
kde_origin_tm_bw_raster <- raster(gridded_kde_origin_tm_bw)
kde_origin_ck_bw_raster <- raster(gridded_kde_origin_ck_bw)
kde_origin_jw_bw_raster <- raster(gridded_kde_origin_jw_bw)8.4.2 Assigning projection system in all raster layer
Code
projection(kde_origin_pg_bw_raster) <- CRS("+init=EPSG:3414")
projection(kde_origin_tm_bw_raster) <- CRS("+init=EPSG:3414")
projection(kde_origin_ck_bw_raster) <- CRS("+init=EPSG:3414")
projection(kde_origin_jw_bw_raster) <- CRS("+init=EPSG:3414")8.4.3 Visualizing the output in tmap
Code
tm_shape(kde_origin_pg_bw_raster) + tm_raster("v") +
tm_layout(legend.position = c("right", "bottom"),frame = FALSE)
Code
tm_shape(kde_origin_tm_bw_raster) + tm_raster("v") +
tm_layout(legend.position = c("right", "bottom"),frame = FALSE)
Code
tm_shape(kde_origin_ck_bw_raster) + tm_raster("v") +
tm_layout(legend.position = c("right", "bottom"),frame = FALSE)
Code
tm_shape(kde_origin_jw_bw_raster) + tm_raster("v") +
tm_layout(legend.position = c("right", "bottom"),frame = FALSE)
9.0 Network Constrained KDE (NetKDE) Analysis Perform NetKDE analysis by using appropriate functions provided in spNetwork package.
9.1 Preparing data 9.1.1 road layer on 4 planning area
Code
pg <- mpsz3414[mpsz3414$PLN_AREA_N == "PUNGGOL",]
tm <- mpsz3414[mpsz3414$PLN_AREA_N == "TAMPINES",]
ck <- mpsz3414[mpsz3414$PLN_AREA_N == "CHOA CHU KANG",]
jw <- mpsz3414[mpsz3414$PLN_AREA_N == "JURONG WEST",]Code
pg_roads <- st_intersection(roads3414,pg)
tm_roads <- st_intersection(roads3414,tm)
ck_roads <- st_intersection(roads3414,ck)
jw_roads <- st_intersection(roads3414,jw)9.1.2 origin data on 4 planning area
Code
origin_pg <- st_intersection(origin_sf,pg)
origin_tm <- st_intersection(origin_sf,tm)
origin_ck <- st_intersection(origin_sf,ck)
origin_jw <- st_intersection(origin_sf,jw)Drop redundant columns
Code
origin_pg <- origin_pg[-c(2:16)]
origin_tm <- origin_tm[-c(2:16)]
origin_ck <- origin_ck[-c(2:16)]
origin_jw <- origin_jw[-c(2:16)]9.1.3 Preparing the lixels objects Before computing NetKDE, the SpatialLines object need to be cut into lixels with a specified minimal distance.
Code
pg_roads = st_cast(pg_roads, "LINESTRING")
tm_roads = st_cast(tm_roads, "LINESTRING")
ck_roads = st_cast(ck_roads, "LINESTRING")
jw_roads = st_cast(jw_roads, "LINESTRING")Code
pg_lixels <- lixelize_lines(pg_roads,
200,
mindist = 100)
tm_lixels <- lixelize_lines(tm_roads,
200,
mindist = 100)
ck_lixels <- lixelize_lines(ck_roads,
200,
mindist = 100)
jw_lixels <- lixelize_lines(jw_roads,
200,
mindist = 100)9.1.4 Gernerating line centre points Next, lines_center() of spNetwork will be used to generate a SpatialPointsDataFrame (i.e. samples) with line centre points
Code
pg_samples <- lines_center(pg_lixels)
tm_samples <- lines_center(tm_lixels)
ck_samples <- lines_center(ck_lixels)
jw_samples <- lines_center(jw_lixels)9.1.5 Peforming NetKDE Puggol
Code
pg_densities <- nkde(pg_roads,
events = origin_pg,
w = rep(1,nrow(origin_pg)),
samples = pg_samples,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 5,
sparse = TRUE,
verbose = FALSE)Tampines
Code
tm_densities <- nkde(tm_roads,
events = origin_tm,
w = rep(1,nrow(origin_tm)),
samples = tm_samples,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 5,
sparse = TRUE,
verbose = FALSE)Choa Chu Kang
Code
ck_densities <- nkde(ck_roads,
events = origin_ck,
w = rep(1,nrow(origin_ck)),
samples = ck_samples,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 5,
sparse = TRUE,
verbose = FALSE)Jurong West
Code
jw_densities <- nkde(jw_roads,
events = origin_jw,
w = rep(1,nrow(origin_jw)),
samples = jw_samples,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 5,
sparse = TRUE,
verbose = FALSE)9.1.6 Visualising NetKDE Insert the computed density values (i.e. densities) into samples and lixels objects as density field.
Code
pg_samples$density <- pg_densities
tm_samples$density <- tm_densities
ck_samples$density <- ck_densities
jw_samples$density <- jw_densities
pg_lixels$density <- pg_densities
tm_lixels$density <- tm_densities
ck_lixels$density <- ck_densities
jw_lixels$density <- jw_densities9.1.7 Rescalling density values
Code
pg_samples$density <- pg_samples$density*1000
tm_samples$density <- tm_samples$density*1000
ck_samples$density <- ck_samples$density*1000
jw_samples$density <- jw_samples$density*1000
pg_lixels$density <- pg_lixels$density*1000
tm_lixels$density <- tm_lixels$density*1000
ck_lixels$density <- ck_lixels$density*1000
jw_lixels$density <- jw_lixels$density*1000Punggol
Code
tmap_mode('view')
tm_shape(pg_lixels)+
tm_lines(col="density", lwd=5)+
tm_shape(origin_pg)+
tm_dots(alpha=0.5)Code
tmap_mode('plot')Observation : According to the Punggol map, it appears that Punggol Way has the highest density. This might be due to the limited convenience of public transportation in the area, leading residents to frequently utilize Grab services for transportation.
Tampines
Code
tmap_mode('view')
tm_shape(tm_lixels)+
tm_lines(col="density", lwd=5)+
tm_shape(origin_tm)+
tm_dots(alpha=0.5)Code
tmap_mode('plot')Observation: According to the Tampines map, Changi Airport stands out as one of the areas with the highest density. This is logical, considering that individuals arriving in Singapore often rely on Grab services, especially since they may have luggage with them.
Choa Chu Kang
Code
tmap_mode('view')
tm_shape(ck_lixels)+
tm_lines(col="density", lwd=5)+
tm_shape(origin_ck)+
tm_dots(alpha=0.5)